
library(xtable)

load(paste("./RDataFiles/072021_btstrap_results_SEED",
            1,".RData",sep=""))
rep.times        = 20
rep.each         = 5
rep.total        = rep.times * rep.each + 1
dim.results      = dim(results); dim.results[2] = rep.total 
dimnames.results = dimnames(results)
dimnames.results[[2]] = paste("btstrap=",1:rep.total,sep="")

results.all = array(NA, dim.results,dimnames.results)
results.all.pmle <- results.all.dr.MLE <- results.all.dr.logit <-results.all

for(i in 1:rep.times){
 
  load(paste("./RDataFiles/072221_btstrap_results_SEED",
             i,".RData",sep=""))
  rep.index = (i-1) * rep.each + 1:rep.each
  results.all[,rep.index,,,] = results[,1:rep.each,,,]
  results.all.pmle[,rep.index,,,] = results.pmle[,1:rep.each,,,]
  results.all.dr.MLE[,rep.index,,,] = results.dr.MLE[,1:rep.each,,,]
  results.all.dr.logit[,rep.index,,,] = results.dr.logit[,1:rep.each,,,]
}
results.all[,rep.total,,,] = results[,rep.each+1,,,]
results.all.pmle[,rep.total,,,] = results.pmle[,rep.each+1,,,]
results.all.dr.MLE[,rep.total,,,] = results.dr.MLE[,rep.each+1,,,]
results.all.dr.logit[,rep.total,,,] = results.dr.logit[,rep.each+1,,,]



results.all.to.Results <- function(results.all){
Results = apply(results.all, 2:3, function(x) x[which.min(x[,"nll time",1]),,])
dim(Results) = c(7,7,rep.total,1)
dimnames(Results) = dimnames.results[c(4:5,2:3)]
Results = aperm(Results,c(3:4,1:2))
return(Results)
}


Results.to.table <- function(Results, subset=F){

b.index = 1:(rep.total-1)

if(subset){
  idx.to.exclude = which(Results[,,"nll time",1] >=0.01)
  b.index = setdiff(1:(rep.total-1),idx.to.exclude)
}

mie = Results[,,"nll time",1]

## Time
mtime = Results[,,"nll time",2]

max.step = "50"

hist(Results[b.index,max.step,"theta","*=0"])
quantile(Results[b.index,max.step,"theta","*=0"], c(0.025, 0.975))
Results[rep.total,max.step,"theta","*=0"]

sd(Results[b.index,max.step,"theta","*=0"])

Results.theta = Results[b.index, max.step, "theta", c(1:2,4:7)]
point.theta   = Results[rep.total, max.step, "theta", c(1:2,4:7)]
test.statistic = t(point.theta) %*% solve(cov(Results.theta)) %*% point.theta
1 - pchisq(test.statistic, length(point.theta))

Results.comp = Results[,max.step,"contrast",1:3]
point.comp   = Results.comp[rep.total,]
sd.comp      = apply(Results.comp, 2, sd)
format(round(cbind(point.comp, point.comp + qnorm(0.025,0,1)*sd.comp, 
                   point.comp + qnorm(0.975,0,1)*sd.comp),3),3)
 
sds = apply(Results[b.index,max.step,,], 2:3, sd)[1:5,]

point.est = Results[rep.total, max.step, 1:5, ]
lower = format(round(point.est + qnorm(0.025,0,1) * sds,2),2)
upper = format(round(point.est + qnorm(0.975,0,1) * sds,2),2)
point.est = format(round(point.est,2),2)

output = t(matrix(paste(point.est, " (", lower, ",", upper, ")", sep=""),5,))
colnames(output) = c("theta","phi","gop","eta 0*","eta 1*")
rownames(output) = c("0","1",".","Household size > 3","Race non-white",
                     "Employed","Married")
return(output)
}


Results.all.dr.MLE <- results.all.to.Results(results.all.dr.MLE)
output.dr.MLE <- Results.to.table(Results.all.dr.MLE,subset = F)
output.dr.MLE.2 <- Results.to.table(Results.all.dr.MLE,subset = T)
# Print the theta
output.dr.MLE[, 1]
output.dr.MLE.2[, 1]


Results.all.dr.logit <- results.all.to.Results(results.all.dr.logit)
output.dr.logit <- Results.to.table(Results.all.dr.logit,subset = F)


Results.all.pmle <- results.all.to.Results(results.all.pmle)
output.pmle <- Results.to.table(Results.all.pmle)
xtable(Results.to.table(Results.all.pmle))

output.pmle[,1]
output.dr.MLE[, 1]
output.dr.logit[,1]
